Noise reduction in chaotic time series by a local projection with nonlinear constraints 
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On the basis of a local-projective (LP) approach we develop a method of noise reduction in time 
series that makes use of nonlinear constraints appearing due to the deterministic character of the 
underlying dynamical system. The Delaunay triangulation approach is used to find the optimal 
nearest neighboring points in time series. The efficiency of our method is comparable to standard 
LP methods but our method is more robust to the input parameter estimation. The approach has 
been successfully applied for separating a signal from noise in the chaotic Henon and Lorenz models 
as well as for noisy experimental data obtained from an electronic Chua circuit. The method works 
properly for a mixture of additive and dynamical noise and can be used for the noise-level detection. 
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I. INTRODUCTION 



It is common that observed data are contaminated by 
noise (for a review of methods of nonlinear time series 
analysis see Q, 0, 0)- The presence of noise can sub- 
stantially affect such system parameters as dimension, 
entropy or Lyapunov exponents ,4]. In fact noise can 
completely obscure or even destroy the fractal structure 
of a chaotic attractor Q and even 2% of noise can make a 
dimension calculation misleading . It follows that both 
from the theoretical as well as from the practical point of 
view it is desirable to reduce the noise level. Thanks to 
noise reduction S S S 113, (IB 113, 113, EE IS E3, E3l 
it is possible e.g. to restore the hidden structure of an 
attractor which is smeared out by noise, as well as to 
improve the quality of predictions. 

Every method of noise reduction assumes that it is pos- 
sible to distinguish between noise and a clean signal on 
the basis of some objective criteria. Conventional meth- 
ods such as linear filters use a power spectrum for this 
purpose. Low pass filters assume that a clean signal has 
some typical low frequency, respectively it is true for high 
pass filters. It follows that these methods are convenient 
for a regular source which generates a periodic or a quasi- 
periodic signal. In the case of chaotic signals linear filters 
cannot be used for noise reduction without a substantial 
disturbance of the clean signal. The reason is the broad- 
band spectrum of chaotic signals. It follows that for 
chaotic systems we make use of another generic feature of 
dissipative motion located on attractors that are smooth 
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submanifolds of an admissible phase space. As results 
corresponding state vectors reconstructed from time de- 
lay variables are limited to geometric objects that can be 
locally linearized. This fact is a common background of 
all local projective (LP) methods of noise reduction. 

Besides the LP approach there are also noise reduction 
methods that approximate an unknown equation of mo- 
tion and use it to find corrections to state vectors. Such 
methods make use of neural networks or a genetic 
programming |l2j | and one has to assume some basis func- 
tions e.g. radial basis functions to reconstruct the 
equation of motion. Another group of methods are mod- 
ified linear filters e.g. the Wiener filter 0], the Kalman 
filter , or methods based on wavelet analysis 0] . Ap- 
plications of these methods are limited to systems with 
large sampling frequencies, and they are confined to the 
neighborhood of every point in phase space. 



The method described in this paper can be considered 
as an extension of LP methods by taking into account 
constraints that occur due to the local linearization of the 
equation of motion of the system. We call our method 
the local projection with nonlinear constraints (LPNC) . 

The paper is organized as follows. In the following 
section we shall present the general background of LP 
methods. The LPNC method is introduced in Sec. IIIII 
and compared with LP methods in Scc. lIVI In Sec.lVlwe 
present methods how to find the nearest neighborhood, 
and examples of noise reduction and estimation are in- 
troduced in Sees . I VII and I VIII In the appendix lA"lone can 
find the multidimensional generalization of the solution 
presented in Sec. IIIII 
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II. LOCAL PROJECTIVE METHODS OF NOISE 
REDUCTION 

Let us consider a scalar time series {x n }, n — 1, 2, N 
corresponding to an experimentally accessible component 
of the system trajectory. We assume that in the presence 
of measurement noise instead of the clean time series x n 
we observe a noisy series x n : x n = x n + ry„ where ry„ is 
the noise variable. The aim of noise reduction methods is 
to estimate the set {i n } from the observed noisy data set 
{x„}, i.e. to find corrections Sx n such that x n +Sx n » x n . 
The corrections 8x n can be estimated on the assumption 
that x n belongs to a clean deterministic trajectory. Let 
us create vectors of the system state x„ using the Takens 
Theorem Q x„ = {x n , x n -r, — , i n -(d-i)r}, where d is 
the embedding dimension, and r is the embedding delay 
that further will be just 1. Now the simple approach is to 
use a linear approximation for the nearest neighborhood 
of a vector x n and then to estimate an unknown 
equation of motion in the embedded space by a linear fit: 
x n+ i = Ax„ + b. The matrix A is the corresponding 
Jacobi matrix and b is a constant vector. 

In LP methods the local linearity of the system dy- 
namics plays the crucial role. The unknown equa- 
tion of motion of a deterministic systems x n +i = 
F (x n , x n —i, . . . , x n -d+i) is equivalent to the presence of 
a constraint H (x n +i> x n , . . . ,x n -d+i) = 0. If the em- 
bedding dimension d is larger than the dimension of the 
attractor then Q constraints appear: 



H q n) (x n+ i, x n ,... ,x n - d +x) = and 



1, 



,Q < d 

w 

where Q depends on the rounded up dimension d a of 
the attractor, Q = d + 1 — d a . Since we apply a linear 
approximation for vectors Xj £ X^^ the constraints 
can be written as 

Hq ^ (Xn-\-l j X n , . . . , X n — d-]-l ) 

d-Q 

a q j ' in) x i - j+1 - q + 6'-W - x i+1 . q = 0, (2) 

where a| and b q ^ n ' > are elements of A and b respec- 
tively. The main problem of LP methods is to find a 
tangent subspace determined by the linear constraints 

Hq 1 ^ (x n +i,x n , . . . , x n -d+i) — and to perform an ap- 
propriate projection on this subspace. Different LP ap- 
proaches make use of different projecting methods, how- 
ever tangent subspaces are found in the same manner by 
all methods, i.e. the subspace should fulfill the condition 

(|2J) and the condition 



~ 2 

Xj, 



[la. [jjj • Since there are several constraints (|2J) and the 
same data will occur in several Takens vectors x ra there 
are many possible corrections 5x n , q to the same observed 
data x n . In the CHS method one makes a compromise 
between different corrections by taking the average 



a 5x r , 

9=1 



where 



Sx r , 



(n) 



(n) 



(3) 



(4) 



is the correction of x„ obtained due to the constraint 
Hq 1 ^ (x n +i,x n , . . . , x n -d+i) = 0. a is some constant < 
a < 1 and \J n H^ n) = X/H { q n) (x n+ i, x n , . . . , x n - d +i) is 
the gradient of the constraint function. 



B. Schreiber-Grassberger method (SG) 

Instead of the perpendicular projection on the sub- 
space defined by J2J) one can perform a projection by 
correcting only one variable If we choose a;„+i_ r as 
the corrected variable where r w d/2 then the corrections 
are 

\**"n-\-l-\-r i ^n+r; ■ ■ ■ j ^n— a-t-l+TV 

Xn — X n (X 

8H S (x n+ i +r , X n + r , . . . , X n -d+l+r )/dx 

(5) 

where s « ^. The approach can be justified as follows. 
If the largest (unstable) Lyapunov exponent is X u > 
and the smallest (stable) Lyapunov exponent is X s < 



we can write 
If A,, 



and 



dx r , 



■■ 1 -Z ^ e |A,|(r-d+l) < 



A s | then the highest precision for determining 
the denominator of the rhs of (JSJl is usually obtained for 
r = d/2: 



mm 

r=0,...,d 
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(6) 



where Ax n is the error connected with the variable x n . 



C. The optimal method of local projection 
(GHKSS) 



A. Cawley-Hsu-Sauer method (CHS) 

The method makes use of a perpendicular projec- 
tion on a subspace corresponding to the constraints @ 



In the GHKSS method [H HJ developed by Grass- 
berger et.al. one looks for a minimization functional 
that fulfills the linear constraints @ by correspond- 
ing corrections received in a one-step procedure. The 
constraints J2J can be written in the equivalent form 
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(a 9:( ™) • y„) + b 9 '(") = where a new vector y„ = 
{x n+ i, x n , . . . , x n ^d+i} is introduced, the dimension of 
which is larger by one than the dimension of the vector 
x„. Vectors a 9 '(") should be linearly independent and 
appropriately normalized, so that multiple corrections of 
the variables are eliminated, i.e. afl^ ■ Pafl = 8 qq i 
where P is the matrix describing the metric of the sys- 
tem. Let Y^^ be a set corresponding to the nearest 
neighborhood of the vector y„. Minimizing the func- 
tional YJ (5xfc • P _1 5xfc for {fc : yfc G Y^^} under the 

k 

above conditions we get a system of coupled equations. 
The next step is to consider all vectors of Y^^ and to 
calculate the average Q 71 ' = | Y jv m | Xk+i, i = 0, 1, . . . , d 

71 k 

as well as corresponding [d + 1) X (d + 1) covariance ma- 
trix 

C *] ] = WNN\ ^2 X k+i x k+ - d" ) £,- n) > ( 7 ) 
1 " 1 k 

here Y^^l means the number of elements in the set. 

Defining Ri = -A= and T^' — RiC^'Rj one can find Q 

orthonormal eigenvectors of the matrix correspond- 
ing to its smallest eigenvalues e q '^ for q = 1, . . . , Q. 

Let the matrix Ilffl = YJ e <h{ n ) e iA n ) define a subspace 

9=1 

spanned by the eigenvectors e q '( n \ Now the corrections 
to the observed signal can be written as follows 

t^iji,: /o(c (8) 

We see that the GHKSS method does not employ mul- 
tiple corrections resulting from constraints (0), but only 
performs a smaller number of corrections following the 
multiple occurrence of the same variable x n in various 
vectors y, : x„ £ y,. 

The solution (JSJ is a generalization of the CHS and SG 
methods. The main difference between the CHS method 
and the GHKSS method is in the subspace of projection. 
While a perpendicular projection of points is used in the 
first case, projection is on a tangent subspace defined by 
the matrix P in the second case. The matrix P should be 
diagonal and such that the first and the last component 
of the vector y n have only small weights e.g. : 

p [0.001 i=0,d, 

I 1 otherwise. 

The efficiency of noise reduction methods can be mea- 
sured by the gain parameter, defined as 

g = m\og(^f^\ (io) 

V °rec2 J 

where cr^ loise = (^(x n — x n ) 2 ^j is the variance of added 
noise and cr^ ed is the variance of noise left after noise 



reduction. The last value is calculated as the square of 
the distance between the vector of noise-reduced data 
and the vector of clean data divided by the dimension of 
these vectors. The definition of the gain presumes the 
knowledge of the clean data X n = {x n }. 

The noise level parameter TV can be defined as the 
ratio of standard noise deviation a no ise to standard data 
deviation a data 

M= <Jnoise_ (n) 
&dat a 

III. THE PRINCIPLE OF LPNC METHOD 

The LP methods described in the previous section 
make use of linear constraints that appear due to lin- 
ear approximation of the system dynamics. Such a linear 
approximation has only a local character and correspond- 
ing coefficients depend, in fact, on the position in phase 
space. If we assume that the nearest neighborhood of ev- 
ery point x„ is characterized by the same coefficients then 
nonlinear constraints appear that can be used for recon- 
struction of the unknown deterministic trajectory. The 
basic advantage of the local projection with nonlinear con- 
straints (LPNC) method introduced here as compared to 
LP methods is its smaller sensitivity to the input param- 
eters estimation. A weak point of the LPNC method is 
its slower convergence rate with respect to the standard 
LP approach. The LPNC algorithm can be accelerated 
but at the cost of decreasing the gain parameter. Like 
other LP methods the LPNC method belongs to the it- 
erative approaches. A single iteration provides only a 
partial noise reduction and a corrected data set serves as 
an input for the next iteration. 

For the one-dimensional case the Jacobi matrix A 
and the additive vector b describing the locally lin- 
earized dynamics at point x n reduce to scalar coeffi- 
cients A = ai (x n ), b = b (x n ), and the linearized equa- 
tion of motion at x n reads x n+ \ = a\ (x n ) x n + b(x n ). 
Let us consider the nearest neighborhood of x n . 

We assume that the set consists of three points 

|i„ , Xk , Xj G X^^ | which are so close to each other that 

their locally linearized dynamics can be approximately 
described by the same pair of coefficients A = a\ (x n ), 
b = b {x n ). When we write down three linear equations 
of motion for x n ,Xk, Xj 

x n +i = ai (x n ) x„ + b (£„) 
Xk+i = ai (x n ) Z k +b (x n ) 
Xj+i = a\ (x n ) xj + b (x n ) (12) 

the coefficients ai(i„) and b(x n ) can be eliminated. After 
elimination we get a constraint that has to be fulfilled by 
the system variables for consistency reasons. 

G [^-nj = x n {Zk+1 ~ Xj+i) + 

x k {Zj+i - i„+i) + Xj (x n+1 - Xk+i) = 0. (13) 
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In the case of a higher dimension d > 1 we have three 
equations of motions but the number of unknown con- 
stants is larger than two, i.e. 

d 

X n +l — 2J a i (Xn) Sn-i+1 + b (x n ) 
i=l 
d 

Xk+l = ^ a i (#n) Sfc-i+1 + b (x n ) 
i=l 
d 

%+i = ^ o» (in) Xj-i+i + b (x n ) , (14) 



where a,i(x n ) are elements of the first row of Jacobi ma- 
trix and b(x n ) is a constant. The corresponding con- 
straint G d for higher dimensional case is as follows 




t-n+l/ 



\t=i 

a i X k-i+lJ ~ %r. 

a^j-j+i^ (5 n +i - ife+i) = (15) 

The extended constraints and the corresponding calcu- 
lations that are valid for all rows of Jacobi matrix A 
are presented in the appendix [S] The condition 113|) 
and l|15fl should be fulfilled for every point x n and its 



nearest neighborhood X.„ N . Similarly as in LP methods 
these constraints are ensured in the LPNC approach by 
application of the method of Lagrange multipliers to an 
appropriate cost function. Since we expect that correc- 
tions to noisy data should be as small as possible, the 
cost function can be assumed to be the sum of squared 
corrections S = J^ s=1 (Sx s ) 2 . 

It follows that we are looking for the minimum of the 
functional 



N 



N 

E 

71 = 1 



X n G c 



(xr) 



(16) 



After finding zero points of 2N partial derivatives one 
gets 2N equations with 2N unknown variables Sx n and 
A„. However, in such a case the derivatives of the func- 
tional (|16|) are nonlinear functions of these variables. For 
simplicity of computing we are interested to pose our 
problem in such a way that linear equations appear which 
can be solved by standard matrix algebra. To under- 
stand the role of nonlincarity let us write the constraint 

G (x^ N ^j in such a way that explicit dependence on the 

unknown variables is seen (the corresponding equations 
for G d (X^ N ) have a similar form) 



G 



(xr) 



— G (X n , X n+ i) + G (5X n , X n+ i) 



G (X^ ViV , SX n+1 ) + G (SX n , 5X n+1 ) . (17) 
Here we introduced the following notation 



G (X n , X„+l) = X n {Xk+l - Xj+l) + Xk (Xj + l - X n +l) + Xj (x n+ l - Xk+l) 

G (<5X„,X n+ i) = Sx n {xk+i - Xj+i) + 5xk (xj+i - x n +i) + 5x 3 (x n+ i - x k +i) 
G (X™,<5X n+ i) ee x n {5x k +i - Sx j+ i) + x k (Sx J+1 - 5x n+1 ) + x 3 (Sx n+ i - 6x k +i) 
G (5Xn,SXn+l) = fan (Sx k+1 - Sx J+ i) + Sx k (Sx j+ i - 6x n+ i) + Sxj (Sx n+ i - Sxk+i) , (18) 



where X^ N = {x„,x k ,Xj}, X n+1 = 
{a; n +i,a:fc+i,a:j+i}, SX n = {Sx n ,Sxk,Sxj}, 
SX n+1 = {5x n+1 ,5x k+1 ,5xj +1 } and x k ,Xj are 
the near neighbors of x n . Indices arc defined as 
{n, j, k : x n ,Xk,Xj £ X^^}. Note that elements of the 
set X n+1 are not necessarily near neighbors to each 
other. 

The approximation we use in l|17|l follows from the 
fact that in general the nearest neighborhood X^ N does 
not include the same indices as the nearest neighborhood 
Xf,i.e. 

{fc:4eXf}^{i:, je Xr}. (19) 

In the case of not correlated noise and under the assump- 
tion that the introduced corrections completely reduce 



the noise effect Sx s = — r/ s (V s= i n) one can neglect 

the nonlinear terms in Eqs. i|18fl i-e. 

G(6X n ,5X n+1 ) = (V„ =1 ,...jv). (20) 

In the equation (|2U[) we use the fact that (rji) — and 

Taking into account the approximation l|2(J|l one can 
write the following linear equation for the problem i|16|) 

M <SX = B, (21) 

where M is a matrix containing constant ele- 
ments, B is a constant vector, and <5X T = 
{Sxi, 6x2, . ■ ■ , fejv, Ai, A2, ■ • • , Xn} is a vector of depen- 
dent variables (T - transposition). In practice it is very 
difficult or even impossible to find the solution of the 
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equation (|21|l for large N. First, it is time consuming to 
solve a linear equation with a matrix 2N x 2N matrix 
for N > 1000. Second, when M becomes singular the 
estimation error of the inverse matrix M ~ 1 is very large. 
Third, we cannot always find the true near neighbors (the 
set X^^) from the noisy data {xi}. Taking into account 
the above reasons it is useful to replace the global mini- 
mization problem (|16|l by N local minimization problems 
related to the nearest neighborhood X^^. The corre- 
sponding local functionals to be minimized are 

S™ = J2 (test + KG d (X™) = mm 

s 

(V„ = i,...,Ar) where 
{s : x s £ or x s G X n+i or ... x s G X„_ d+ i}(22) 

We can consider the minimization problem H22|) as a cer- 
tain approximation of Functionals H22|) are linked 
to each other due to the fact that the same variable Sx n 



appears in 6-d different minimization problems (|22|l . The 
global problem 116JI is equivalent to Eq. Ij21|l with 2N un- 
known variables that should be found single-time. The 
problem (|22|l is equivalent to a system of coupled equa- 
tions that should be solved several times and as a re- 
sult one gets an approximate global solution. Writing 
Eq. (|22|l in the linear form i.e. calculating zero sites of 
corresponding derivatives and using Eq. I|20|) one gets N 
linear equations as follows 

M„-5X*=B n (V n= i,... lJV ), (23) 

where (SX*) T = {Sx n , Sx k , Sxj, 6x n+1 , Sx k +i, Sx j+1 , A„}. 
The matrices M„ corresponding to (|22|l avoid the dis- 
advantages of Ij21|) . i.e. they are not singular, their 
dimension is smaller and they do not substantially 
depend on the initial approximation of near neighbors. 
The matrix M„ for one-dimensional case is given by 





2 

















Xk+1 


- X j+ i 
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X j+ i 


— X n +\ 
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X n +1 


- x k +i 


M„ = 











2 








X J 


- x k 
















2 





J a 


— Xj 
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Xk 


Xn 




. X k +1 - Xj + \ 




-1 x n +\ — X k +1 


Xj ■ X~fa 












(24) 



Vector B„ has the form = 
{0,0,0,0,0,0,-G(X^,X n+1 )}. 

IV. COMPARING LPNC METHOD TO LOCAL 
PROJECTION METHODS 

Let us illustrate the LPNC method by taking into 
account the cost functional (|22|l (it will be written as 
S LPNC ) 

s lpnc = J- 5x1 + A[«n (Xk+1 - x 3 +i) + 

i 

x k (xj+i - x n +i) + xj (Z n +i - x k +i)] = mm . (25) 

The corresponding cost function S GHKSS that is used in 
the standard local projection method e.g. in the GHKSS 
method [3 is 

s ghkss = 5x} + Ai (x n a + b - x n +i) + 

i 

A2 (xjd + b — Xj+i) + A3 {Xk<x + b — Xk+i) = min . (26) 

If we were in the position to find exact solutions for the 
minimization problems l|25|l and (|26|l then both results 



would be the same since (|2*3|) can be obtained from lt2T>|) 
after elimination of the parameters a and b. 

In both cases the variables ^x n ,Xk,£j G X^^j be- 
long to the nearest neighborhood of the variable x n . The 
index i = ^k,k+p : Xk G X^^l runs through all in- 
dices of the variables appearing in (|25() and (|26|l while 
the variable Xk+ P = Ixi ■ %i-p G X^^j corresponds to 

the p-iterate of Xk ■ Parameters a and b can be calculated 
from a linearized form of the equations of motion at the 
point x n . 

In practice the minimization problems S LPNC and 
gGHKSS are no ^. e q U i va i en t because in both cases dif- 
ferent approximations are used. These differences are: 
(i) Eq. (|25|l is nonlinear against corrections 5xi. In this 
case the approximation consists in a linearization, (ii) 
For Eq. I|26l) the exact values of the parameters a and 
b are unknown. The approximation means that a and b 
are estimated from noisy data. 

Fig. Hand [21 present a comparison between results re- 
ceived by the GHKSS and LPNC methods. Fig.[T]shows 
that the gain parameter Q depends on the number of 
neighbors, which is an input parameter of both methods. 
One can see that for LPNC method the gain parameter is 
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FIG. 1: The plot of the gain parameter Q versus number of 
iterations of the GHKSS method (squares) and LPNC method 
(triangles). Lorenz system TV = 78%, N = 1000. 
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FIG. 2: The plot of the gain parameter Q versus number of 
neighbors of the GHKSS method (squares) and LPNC method 
(triangles). Lorenz system TV = 78%, N = 1000. 
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FIG. 3: The random data from uniform distribution 
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FIG. 4: The random data shown in the Fig. after noise 
reduction with the LPNC method 



THE NEAREST NEIGHBORHOOD 
ASSESSMENT 



more robust to changes of the number of neighbors than 
for the GHKSS method. In Fig.[3the dependence of the 
gain parameter on the number of iteration steps of the 
methods is shown. One can see that LPNC method fin- 
ished reduction at the maximal efficiency what is not the 
case of GHKSS method, so the former method is easier 
to use since it does not need estimation of the iteration 
number. 

If we consider uniformly distributed stochastic vari- 
ables (see Fig. |3J) the LPNC method reduces the noise 
very well, and as a result all data are represented as a 
neighborhood of a point attractor (see Fig. 0} while a 
complete noise reduction would correspond to a phase 
portrait consisting of a single point. In fact, for the case 
considered we observed for the LPNC method a noise 
reduction of about 96% of data variance. 



The LP methods are local. It follows that features of 
the nearest neighborhood of every point x n in the 

phase space play an important role. Usually the near- 
est neighborhood is estimated by the smallest distance 
approach that makes use of the standard Euclidian ge- 
ometry. We have found, however, that our LPNC method 
works much better when the Delaunay triangulation ap- 
proach [2l| is applied for the nearest neighborhood esti- 
mation. 



A. The smallest distance approach (SD) 

In the smallest distance approach the Euclidian met- 
ric is used, i.e. first the distance between every pair of 
points in the Takens embedded space is calculated as 



\J {Xi - Xjf + . . . + (xj_ (d _i) T - Xj_ (d _i) T ) 2 and 
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FIG. 5: Delaunay triangulation for a set of nine points in 
a two-dimensional space. The near neighbors are connected 
by bold lines. Sets T (x;), i = 1,2, . . . , 9 are limited by thin 
lines. 



then the nearest neighborhood of a point x„ is de- 

fined as a set of v points fulfilling the relation 



{ Xj exr,V fe x^X^ 



NN 



(27) 



Let us stress that this definition depends on the chosen 
value of the v parameter, i.e. on the assumed number of 
near neighbors, v = 2, 3, 



B. The Delaunay triangulation approach (DT) 

To find the nearest neighborhood relations for the 
LPNC method we have used the Delaunay triangulation 
pl|. In general the triangulation of any set of points 
X = {x,;} e l d is a collection of d-dimensional sim- 
plices with disjoint interiors and vertices chosen from X. 
There are many triangulation of the same set of points 
X. One of the best known is the Delaunay triangulation 
(see Fig. EJ. Let T(x n ) be a part of the space M. d that 
contains all points that are closer to x n than any other 
point Xj from the set X 

T(x„) = {zel j ,x, eX : V^nllz-XnH < \\z-XjW}. 

(28) 

If rriij = (xj +Xj) /2 belongs to both sets T(xj) and 
T (x.j ) then by definition the point Xj is the nearest neigh- 
bor of Xj received due to the Delaunay triangulation. By 
the above definition every point x„ belongs to its nearest 
neighborhood x„ £ X^^ . 

In practice the Delaunay approach can be performed 
as follows. A pair of points x„ and Xj are near neighbors 
provided that there are no other points X& (k ^ j, n) 
belonging to the hypersphere centered at the point m ra j 
and of the radius r n j = ||x„ — Xj || /2. 

In Fig. [S] two cases are presented when in the two- 
dimensional space a) the point Xj is not the nearest 
neighbor of x„ and b) the point Xj is the nearest neighbor 
of x„ . 




b) 



FIG. 6: Illustration of nearest neighborhood search by DT 
approach a) Xj and x n are not near neighbors, b) Xj and x n 
are near neighbors 



T 



X 




FIG. 7: Chaotic Henon map without noise 



The DT method has the advantage that triangles ap- 
pearing due to connections of near neighbors are almost 
equiangular (see Fig.^J). This property is the main reason 
for using the DT method in search of the near neighbors. 
The disadvantage of this method is a slowing down of 
numerical calculations. 



VI. EXAMPLES OF NOISE REDUCTIONS 

The LPNC method has been applied to three systems: 
the Henon map, the Lorenz model |22( and the Chua cir- 
cuit I2II2II23. Figures El present the chaotic Henon 
map in the absence and in the presence of measurement 
noise as well as a result of the noise reduction. Table [I] 
presents the values of the gain parameter for the Henon 
map and for the Lorenz system. 

To verify our method in a real experiment we have 
performed the analysis of data generated by a nonlin- 
ear electronic circuit. The Chua circuit in the chaotic 
regime [23l [24| has been used and we have added a mea- 
surement noise to the outcoming signal. The noise (white 
and Gaussian) came from an electronic noise generator. 
Figures - ^] show a clean signal coming from this 
circuit, the signal generated by Chua circuit with mea- 
surement noise (M — 96.5%) and the same signal after 
the noise reduction with the LPNC method (G = 6.38). 
Table ITT1 presents values of the Q parameter and the per- 



-4-3-2-101234 
X 

n 

FIG. 8: Chaotic Henon map with a measurement noise Af : 
69%. Note the difference in scale. 



1.9 2.0 2.1 2.2 2.3 2.4 

FIG. 10: The stroboscopic map corresponding to a clean 
trajectory in the Chua circuit. 




FIG. 9: Chaotic Henon map with a measurement noise after 
noise reduction with the LPNC method, Q — 4.3 



centage of eliminated noise for several values of the noise 
level in the Chua circuit. 

All the above applications of the LPNC method con- 
sider the case of measurement noise that has been added 
to the signal in numerical or electronic experiments. 
However, our LPNC method can also be applied to dy- 
namical noise i.e to the noise which in experiments is in- 
cluded in the equations of motion |2rjj | . In such a case one 



TABLE I: Results of noise reduction by the LPNC method 
for the Henon map and Lorenz model. 



System 



Q percent, of eliminated noise 



Henon N = 1000 10% 9.58 

Henon N = 3000 10% 10.04 

Henon N = 1000 66% 5.08 

Lorenz N = 1000 78% 5.85 

Lorenz N = 3000 76% 6.02 

Lorenz N = 1000 34% 7.21 



89% 
91% 
69% 
74% 
75% 
81% 




FIG. 11: The stroboscopic map received from the Chua cir- 
cuit in the presence of a measurement noise N = 96.5%. Note 
the difference in scale. 




+ 2.1 



1.9 2.0 2.1 2.2 2.3 2.4 

FIG. 12: The stroboscopic map received after the noise re- 
duction by the LPNC method applied to data presented at 
Fig, till 
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TABLE II: Results of noise reduction by the LPNC method 
for the Chua circuit with a measurement noise (N — 3000). 



N 



percent, of eliminated noise 



24.9% 
28.3% 
46.1% 
73.7% 
90.6% 
96.5% 



5.4 
4.9 
7.0 
4.81 
7.4 
6.4 



71% 

68% 
80% 
67% 
82% 
77% 




2.1 2.2 2.3 
V [V] 



FIG. 13: Stroboscopic map received from the Chua circuit in 
the presence of a mixture of a measurement and dynamical 
noise N w 22%. 



cannot compare the noisy data with the clean trajectory 
since the latter one does not exist anymore, and there 
are only e-shadowed trajectories § that can be approxi- 
mated by means of the LPNC method. FigurelTSIshows a 
measured signal generated by a Chua circuit where a mix- 
ture of measurement noise and dynamical noise occurs. 
Figure IbH shows the result of noise reduction applied to 
such a signal. 



VII. NOISE LEVEL ESTIMATION BY LPNC 
METHOD 




FIG. 14: Stroboscopic map received after noise reduction by 
the LPNC method applied to data presented at Fig. ED 



TABLE III: Noise level estimated by the LPNC method for 
the Chua circuit with measurement noise (N = 3000). 



N 


CTnoise [mV] 


O noise [iH 


0% 





5.5 


3.1% 


30.4 


28.9 


6.2% 


60.8 


53.7 


12.3% 


121.7 


110 


24.9% 


243.4 


235 


28.3% 


304 


305 


46.1% 


486 


454 


73.7% 


973 


938 


90.6% 


1520 


1375 


96.5% 


2120 


1844 



V] 



and the fact that the method can be used only for low- 
dimensional systems. On the other hand the LPNC 
method can be applied for estimation of any noise level 
including a large one. In Table UTTI we have presented the 
estimated noise level a no i se for the Chua circuit. 



The LPNC method introduced in the previous section 
can be used to quantify the noise level of data. The noise 
level, i.e. the standard deviation in noisy time series, may 
be approximated as the Euclidian distance between the 
vectors {xi} and {xi} representing the time series before 
and after noise reduction 0] 



\ 



1 N 

N 



Xi) 



(29) 



i=l 



The main disadvantage of the LPNC method used for 
the noise level estimation is its small rate of conver- 
gence with respect to other known methods pl. l27l.l28l . l29]] 



VIII. CONCLUSIONS 

In conclusion we have developed a method of noise 
reduction that makes use of nonlinear constraints which 
occur in a natural way due to the linearization of a deter- 
ministic system trajectory in the nearest neighborhood of 
every point in the phase space. This neighborhood has 
been determined by Delaunay triangulation. The method 
has been applied to data from the Henon map, Lorenz 
model and electronic Chua circuit contaminated by mea- 
surement (additive) noise. The efficiency of our method 
is comparable to that of standard LP methods but it is 
more robust to input parameter adjustment. 
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APPENDIX A: MULTIDIMENSIONAL VERSION 
OF LPNC METHOD 

In Sec. IIIII the LPNC method has been presented for 
one-dimensional systems. Here we show the generaliza- 
tion of this approach for d-dimensional dynamics. For 
one-dimensional problems the Jacobi matrix of the sys- 
tem does not appear explicitly in our method. For higher 
dimensional models the corresponding Jacobian A has to 
be calculated but we manage to minimalize errors occur- 
ring by its estimation. The linearized equation of motion 
for vectors from the nearest neighborhood X^ of a vec- 
tor x„ can be written in the form 

x„ + i = A • x„ + b. (Al) 

In such a case one needs three vectors x„, x/., x, € X^ w 
to write constraints corresponding to Eq. H13|) . In com- 
parison with to the one-dimensional case the number of 



near neighbors i.e. the number of points in the set X„ 
must be larger to allow a unique estimation of the Jaco- 
bian A. We assume that the Jacobi matrix can be ap- 
proximately received by minimalization of the following 
cost functional 

^C^ml^s O m 2X s — 7- -(-... 
s 

• • • + aw^s-(d-i)T ~ ^s+i-mr) 2 = rain 

(V m= i,..., d ) and { S :x s eX™}, (A2) 

where = [A]y. By analogy with Eq. (|13|) we introduce 
the follow 



(^-nj = a mlG (lt^ N , X„+l_ mr j + 
a m2G ^X„_ r , X„+i_ mT ^ + . . . 

G (x n _ (d _ 1)T ,X„ +1 _ mT ) = 

(V m= i,..., d ), (A3) 

where we used the notation corresponding to equa- 
tion (H2J i.e. 



G (X% N , X n+ i) = x n (x k +i - x j+ i) + x k (x j+ i - x n+ i) + xj (x n+ i - x k +i) 
G (X„_ s , X n+t ) = x n - s (x k+t - Xj + i) + x k _ s (x j+ i - x n+l ) + Xj- S (x n+ i - x k+l ) , (A4) 



X n+S = {x n+s ,Xfe +s ,x J+;j } (V s =o,±i,±2,...) where 
|n, j, k : x„ , Xfe , Xj £ 

x « N j- Since the clean trajectory 
is not known thus in the Eq. (|A3|) the observed variables 
X^ N , X„ + i_ mr etc. are used. 

In such a way the equation l|20|l can be written in a 
more general way as 

d 

a mlG (<5X„_ <5X„ +1 _ mT ) = 

2 = 1 

(V n= i,...,jv,V m=li ... )rf ) (A5) 

where we use 

G (SX n - s ,SX n+ i) = 5x n - s (5x k +i - Sxj+i) 
+5x k - s (Sxj + i - 6x n+ i) + 5xj~ s [5x n +i - 5x k +i) (A6) 

SX n+s = {Sx n+s ,Sxk+ s ,Sxj +s } (V s=0 ,±i,±2,...), 
where {n, j, fc : x„, x fe , x 3 <E X^}, 8x n = 
{Sx n ,5x n -T, ...,fo„_( d _i) r } and {n:x„eXj r!f }. 



Now the cost problem (|22H can be transformed to the 
form 

= E (<^) 2 + X nG d m (X ? r ) = min 

s 

\^n=l,...,N:^m=l,...,d) 

and {s : x s G X% N or x s 6 X„ +1 } . (A7) 

Finding zeros of partial derivatives of the functional (|A7|) 
one can linearize this problem and write it in the form 
similar to the Eq. $Z3\i 

M„-(5X J A l = B„ (V n=1 ,... iJV ). (A8) 

Vectors 5X^ and B n occurring in Eq. (|A8|) are equal to 

( SX n) T = {^n-(d-l)T.fen-(d-2)r)--- 

. . . , 6x n , 5x n+ x, A* , Xl,..., A^}, (A9) 

Bl = {0 I 1 ... 1 0,-Gf(xr),-^(Xf) J ... 

...,-^(X^)XA10) 
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where the number of zeros do appearing in depends 
on the values of r and d and for the case r = 1, do = d+1. 



J 



Elements of the matrix M„ can be written as 



l M n] mm 

[M«] lm 

[ M n]lm 



2 (V m= i v .. i d ) 

d 



[M„] m , 



/^am8 (xk+i - Xj+i) (V m =d +i,...,d 0+ d+i) and {7 : x/ G x„} 

s=l 
d 

i+i — x, i+ i) (V m =da+i,...,(io+d+l) an d {Z : X/ G x^} 

s=l 
d 

^ ^ ^ms (Xn+1 ) (y m =d +i,...M +d+i) and {Z : x; G Xj} 

s=l 
d 

) (V m =d Q +i,... )( i +d+i) and {Z : x/ G x„ + i} 

s=l 
d 

^a ms (x n _ mT — Xj- mT ) (V m= d D +i ) ...,d +rf+i) and {Z : x/ G x^+i} 
s=l 
d 

) (V m= d +i i ... ! d n +d+i) and {I : x; G Xj + i} 



where the remaining M„ elements vanish and x/ G x„ 
means that the variable x; is a component of the x„ vec- 
tor. 

The operator + = in lAllf) has the same meaning as 
in the programming language C++, i.e. if elements of 



(All) 



r 



the matrix [M„] m/ occur in a few places (e.g. : x n G 
x„ and x„ G x„ + i Vd>i, V T= i) then the elements at 
the rhs of such equations have to be summed up. 
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